Electronic damping of molecular motion at metal surfaces 
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A method for the calculation of the damping rate due to electron-hole pair excitation for atomic 
and molecular motion at metal surfaces is presented. The theoretical basis is provided by Time 
Dependent Density Functional Theory (TDDFT) in the quasi-static limit and calculations are per- 
formed within a standard plane-wave, pseudopotential framework. The artificial periodicity intro- 
duced by using a super-cell geometry is removed to derive results for the motion of an isolated 
atom or molecule, rather than for the coherent motion of an ordered over-layer. The algorithm is 
implemented in parallel, distributed across both k and g space, and in a form compatible with the 
CASTEP code. Test results for the damping of the motion of hydrogen atoms above the Cu(lll) 
f) |. surface are presented. 
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Considerable progress has been made in recent years in understanding the fundamental processes involved in gas- 
surface interactions. This has been based on the parallel developments of large-scale electronic structure calculations 
based on density functional theory, combined with multi-dimensional quantum and classical analysis of the dynamics 
Despite these advances there remains one key area that is still largely unexplored and poorly understood; the 
process of energy dissipation into substrate degrees of freedom. Although this is known to be of central importance 
q ■ in many situations |2|, there exist no 'real' calculations to date for the energy loss to either phonons or electrons in 
" ^ | the surface. 

In particular there have been a number of recent experiments that have shown convincing evidence that energy 
dissipation by the creation of electron-hole pairs is a significant effect in gas-surface dynamics. Gostein et al [3| 
Qh carried out a detailed state-to-state analysis of H2 scattering from Pd(lll) and showed that, for example, in the 
vibrational relaxation of (y = 1, J = 1) to (y = 0, J = 5) an average of 120 meV is lost to the substrate during 
the scattering event, presumably to electron-hole pair formation. Nienhaus and co-workers Q measured directly the 
hot electrons and holes created at Ag and Cu surfaces by the adsorption of thermal hydrogen and deuterium in the 
form of 'chemicurrents' in a Schottky diode. Finally, Huang et al [a] have studied NO scattering from Au(lll) and 
have concluded that the main sink of energy for the vibrational relaxation of v = 2 molecules is the surface, with the 
strong dependence of the de-excitation probability on incident energy providing evidence that an electron-hole pair 
mechanism is the dominant factor. 
0^ ' We carry out a calculation of the ground state properties of an interacting surface/molecule system using a plane- 
wave basis and a super-cell geometry, and use these results to evaluate the friction coefficient associated with the 
motion of a molecule at a chosen position and in a direction of choice. This is achieved using the well established 
'Golden Rule' expression [|| 0, H[ that may be obtained by applying Time Dependent Density Functional Theory 
(TDDFT) together with a quasi-static limit [3], or less stringently by applying the Golden Rule directly to the 
k> I available Kohn-Sham states [y]. Essentially the theory is as described by Hellsing and Persson Q or Liebsch Q. The 
super-cell method has the advantage that it retains the continuous spectrum of one-electron excitations, unlike cluster 
models and this is important for the interactions considered here. In first principles calculations of molecule- 
surface systems a super-cell of sufficient size is usually chosen to prevent any significant interaction between the 
adsorbates in neighbouring super-cells. When considering electron- hole pair excitation a slightly more subtle effect 
must be taken into account, arising from the enforced periodicity of the perturbation that produces the electron-hole 
pairs. We are primarily interested in the energy loss by electron-hole pair excitation due to the motion of an isolated 
molecule interacting with the surface, whereas the super-cell geometry will naturally describe the damping of an 
ordered over-layer. For the periodic system, conservation of crystal momentum prevents transitions occurring that 
will occur for an isolated molecule interacting with the surface. Results for an isolated molecule are derived from the 
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available periodic perturbation, and a significant difference is found between the energy loss behaviour of the periodic 
and isolated systems. 

To test our method we investigate the friction coefficient of an H atom above the hep hollow site of a Cu(lll) surface. 
Spin is included explicitly in the Kohn-Sham theory using the gradient corrected local spin density approximation for 
exchange-correlation (LSDA-GC). In the next section the evaluation of the dynamic self-energy and friction coefficient 
from Kohn-Sham results using a plane- wave basis is described. In section HTT1 the implementation of this as a parallel 
algorithm is described, along with a brief description of the performance of the algorithm. Results for H/Cu(lll) are 
discussed in section HVl and section [V] is the conclusion. 



II. THEORY 



The experimentally measurable energy loss spectrum for a particular mode is directly related to the dynamic self- 
energy, A(u>). Since we are interested in the energy loss, the imaginary part of this self-energy is the required quantity 
and this can be expressed using Q 

Im AH = ^EE |(^l^ ff (r,c)|<''; 

k,k' n,n' 

x(f(^)-f(4))s(^-eZ + 4) (1) 

where ipk an d e k are the Kohn-Sham wavefunctions and energies resulting from a density-functional description of 
the ground state, and /(e) is the Fermi-Dirac occupation function. Spin degeneracy is assumed in Eq. |T]), hence the 
factor of 2; the extension to spin polarised systems is straightforward. The effective field, <jf s , is the TDDFT effective 
field with contributions from the field of the displaced nuclei, the Coulomb field of the induced charge density and 
a contribution from exchange-correlation. Equation fl} can be derived by direct application of the Golden-Rule and 
a single electron approximation Q , or through TDDFT with the assumption of a slow time-dependent perturbation 
(see H and [lfj for more clarification of the role of TDDFT). 

Equation (JTJ) provides the rate of energy loss of a mode of frequency oj due to the excitation of electron-hole pairs 
as t _1 = Im A(u>)/u>. Taking the quasi-static limit oj — > results in the rate of energy loss for the motion of the atom, 
and this can be expressed as the friction coefficient rj defined via the standard Langevin equation. For motion in the 
direction h, r) is given by (see Hellsing and Persson Q, or Plihal and Langreth (Tlj ) 

77 = M lim^^o I m A(w) jijj 

= 27rfiEk,k'E„,„'|(^|h.VV|^:)| 2 ''S(e F - e£)5(e F - 4) (2) 

where e-p is the Fermi energy, M is the mass of the molecule, and h.W is the static limit of —etjf given by the 
derivative of the Kohn-Sham self-consistent potential in the direction h. 
Equation ^ is evaluated using wavefunctions of the form 

<(r) = u£e ik ' r = ^=E Q (k)e i(k+g) ' r (3) 

where Vq is the volume of the super-cell, and only the wavefunctions sampled at a finite number of k points in the 
irreducible wedge of the Brillouin Zone are available. The change in the Kohn-Sham potential due to the motion of 
an isolated atom is not directly available from a super-cell calculation, but we can obtain the change due the motion 
of a periodic lattice of atoms, SV^ icc , as a finite difference. 



A. Obtaining the interaction of an isolated molecule from that of an over-layer 

To calculate the dynamic interaction of an isolated molecule with a surface within a super-cell geometry care must 
be taken with the interpretation of Eq. (JTJ) . If h.W is obtained directly as a periodic function, drastic consequences 
result - the integral becomes zero for k ^ k'. This is due to the super-cell geometry enforcing a conservation of crystal 
momentum that would be correct for describing the interaction of a real over-layer of molecules in coherent motion, 
but is not physically realistic for the aperiodic single molecule/surface system that concerns us here. In light of this 
we must obtain the change in the Kohn-Sham potential due to the motion of an isolated molecule, c> l^|^f ated , from 
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the change due to the coherent motion of an over-layer, 8V^ icc . The relationship between the two can be found by 
applying linear response theory Although the theory is given here for a local pseudopotential, the generalisation 
to non-local pseudopotentials is straightforward. 

We begin by considering the change in the total pseudopotential, S<p(r), due to an infinitesimal change in the 
position of the atom in each super-cell, (5u x (l). With the pseudopotential due to an atom at x denoted V P seudo(r — x) 
this takes the form 

S(j)(r) = W P scudo(r - x - l).Su x (l), (4) 
1 

where 1 is a lattice vector. The change in the Kohn-Sham potential can be obtained from S<p via the static inverse 
dielectric function, %g(r, r') that corresponds to the original Kohn-Sham calculation [13j . This gives 

5V£L*(r)= J £K i(r,r')^(r')d 3 r' (5) 

or 

^ ice (r)=^R x (r,l)Ju x (l) (6) 
l 

where 

Rx(r, 1) = J e K l{r, r')VV pscudo (r' - x - l)d 3 r' (7) 

is the derivative of the Kohn-Sham potential with respect to the change in the position of an atom at x + 1. 

For an isolated atom we require <5V^ atcd (r), the change in the Kohn-Sham potential due to the change in the 
position of an atom at x. The formally correct way to obtain this is to obtain e^g(r, r') from the original DFT 
calculation (for an example of this type of calculation see Godby et al [H|]), obtain R x (r, 1), and use this to evaluate 

^Xted(r) = Rx(r,0).5u x (0). (8) 

However, provided R x (r, 1) is well localised around the site of the perturbed atom such that there is little overlap 
between the responses to the motion of atoms in adjacent unit cells we may take 

^ LcdW « 6(r)^R x (r,l).<5u x (l) (9) 
l 

« e(r)SV£ t | ce (r) 

where 0(r) = 1 within a Wigner-Seitz cell centred on the perturbed atom, and is zero elsewhere. 

From the super-cell calculation the Kohn-Sham potential for the atom at x ± h is evaluated to obtain h. W as 



h.W = 



sv. KS 

isolated 



5|u x (0)| 

^ ([^lattice - e F]x+h - [lattice ~ e F]x-h) (10) 



2|h 

where the variation of the zero of the Kohn-Sham potential has been corrected for using the Fermi energy associated 
with the self-consistent results for each atomic position. 

Equation (|10[) will be accurate providing the response is isolated within the Wigner-Seitz cell centred on the 
perturbed atom (or, equivalently, Vj^.f icc (r) is negligible at the border of this Wigner-Seitz cell), and |h| is small 
enough. In practice Eq. (|10p corresponds to reducing the volume of integration in Eq. ([1]) from the entire lattice to 
one Wigner-Seitz cell. In section lTvl the consequences of considering the motion of an isolated atom, as described above, 
are investigated by comparing results with those obtained by treating the motion as that of an ordered over-layer of 
atoms. 



B. Evaluation for a plane- wave basis 



A plane-wave calculation results in a set of Kohn-Sham states on a finite grid (in real space or reciprocal space) 
sampled at a finite number of k points and for a finite number of bands, hence we must obtain the discrete equivalent 
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of Eq. ([2]) to evaluate 77. It is also desirable to apply the space-group symmetry to reduce the number of k points 
that must be considered. This is achieved as follows. 

First Eq. fl} is discretised using a uniform grid of k points throughout the 1 st Brillouin Zone, and a conventional 
smearing function 

Im AH = E \mr S M\^) 



k,k' n,n' 



(11) 



where N is the number of k points, /(e) is the Fermi-Dirac occupation function, and —S(x) is the derivative of the 
'squashed Fermi-Dirac' 15] function used for occupation numbers in the original Kohn-Sham calculation. 8(x) is 
given by 



xr \ 1 y ( l 2 



" At 



(12) 



where A is an energy parameter describing the width of the function, and t ■■ 
To obtain r\ we take lim w _»o I m resulting in 



2 ^EE \m^y\i>i) 2 (m) - f(4i) 



k,k' n,n' 



(13) 



where —S'(x) is the 2 nd derivative of the 'squashed Fermi-Dirac' function, and we have taken the static limit of (jf s . 
Equation (|13p reduces to Eq. ([2} in the zero temperature and continuum limit. It is important to note that Eq. (|2|) 
cannot be directly discretised by replacing the 8{x) functions with S(x) as this corresponds to taking the continuum 
and zero temperature limits before discretisation. This would result in the inclusion of a contribution to r\ from 
transitions that should not be present (such as from a state to itself). 

Next we reduce the k points that must be considered to those within the irreducible wedge of the Brillouin Zone 
using the space-group of the lattice. Since the set of k points possess point group symmetry Eq. (|13p immediately 
takes the form 



1 IW 2 

= 2nh ^ E E E | <v^ k |h.wi< k ,)[ (/c 



s S,S'k,k'n,n 



(14) 



where IW denotes a sum over points in the irreducible wedge, and S,S' are space-group operators composed of a 
unitary transformation, P, and a non-symmorphic translation, w. The number of operators in the space-group is 
denoted by N$, and — N-^/N where is the number of distinct k points in the entire Brillouin Zone created 
by applying the complete set of Ns space group operators to the point k in the irreducible wedge. There are no 
symmetry operators associated with the eigenvalues since eJJ = £p k . Transforming the integrand in real space by S" -1 
and using the identity [l6l ] 



results in 



«(r)=^(Pr + w)=^) 



, IW 2 

'? = 2 * h W s EEE |W$k|h.W|V#>| (/(e£) ~ f(4)) 



(15) 



S k,k' n,n' 



x<5'(e£ - e£,)u; k -u; k / 



(16) 



where repeated sums over the same operator are removed and h. W is required to possess the space-group symmetry 
corresponding to the available k points. This requirement means that the original Kohn-Sham calculations must be 
performed with symmetry low enough to allow the motion of the atom in the direction we are interested in, even 
if the instantaneous position of the atom corresponds to a higher symmetry. In practise we also only calculate the 

quantities within the sum where the f/(e k ) — /( e k')) ^'( e k — e k') factor is greater than some small value, preventing 



FIG. 1: For an atom at O integration must be performed over the Wigner-Seitz unit cell centred at O (dashed line). For the 
parallelogram (which represents the originally chosen unit cell) the function h.VV is not localised within the unit cell, and 
does not possess space-group symmetry. 



the cost of calculating matrix elements that make negligible contribution to r\ and limiting the sum to a finite number 
of bands. 

For the the motion of an over-layer of atoms, where h.VV is periodic, a similar treatment results in 



IW 2 

r? k=k < = E |<^|h-W|V^'}| (JK) - /(£')) S'& ~ e£V (17) 

k 71,71' 

where only k = k' transitions contribute. 

In order to obtain the matrix elements we require the Kohn-Sham states corresponding to the images of the 
available k points under transformation by a point group operator, i/'pk- Transforming ip£ in reciprocal space gives 
us the coefficients of each ipp k , and this is achieved by applying Eq. (p~5|) to Eq. ([3]) and projecting out the required 
coefficients to give 

C"(Pk) = C£_ lg (k)e- i ( Fk+ «>- w . (18) 

The wavefunctions in real space are then obtained using a FFT 

ul = FFT[Cg (k)] (19) 

and the required integrals are evaluated as sums over a unit cell (this is analytically correct for a plane- wave basis). 
This cell must be chosen such that h..W is negligible at its borders and such that the 'truncated' potential possesses 
the correct space-group symmetry. As mentioned previously the ideal choice is a Wigner-Seitz cell centred at the 
interacting atom, as is shown schematically in Fig. Q] for a 2d hexagonal unit cell with the atom of interest at O. 
To perform the sum it is necessary to divide the integration volume into different regions Xj such that transforming 
region Xi by a vector constructs the Wigner-Seitz cell from the original unit cell. The matrix element then takes 
the form 

(^ k |h.W|^;) = E^E^fe)^)^^)^-^'!'^) (20) 

where h.VF and the wavefunctions are available over a grid of N r real space points, rj. 

It should be noted that for an isolated atom the integral is carried out over one unit cell as a consequence of h.VV^ 
being localised in that cell, whereas for the coherent motion of an over-layer the integral is carried out over one unit 
cell as a consequence of the periodicity of the integrand. Hence, to evaluate ??k=k' the rearrangement of Eq. (|20[) is 
not necessary. 



6 



{g}i 



{g} 2 



{k}i 



1 



l + N G 



{k} 2 



2 



2 + N G 



N G x P{g} 



FIG. 2: Distribution of data across processors - each box is one processor and contains the processor number. Columns are 
processors that hold data for subregions of real or reciprocal space, and rows are processors that hold data for subsets of k 
points. 



We are interested in large scale systems, hence a parallel implementation of both the original Kohn-Sham calculation 
and the evaluation of r\ is desirable. Multiprocessor algorithms for plane-wave Kohn-Sham methods are available, so 
we describe the latter only. 

First we consider the distribution of data between the N p available processors. These are divided into Ng groups 
of -P{g} processors, and the k points are distributed equally between these groups of processors. Data specific to a k 
point is stored on the associated group. The g points are then distributed across the P{g} processors in each group, 
and data associated with each g point held on the associated processor. Data that is not dependent on k is the same 
for each group but distributed across processors within the group, and data that is not dependent on g is the same 
within each group but distributed across groups. This distribution is shown figuratively in Fig. [2| where the k points 
in the n th group are denoted {k}„, and the g vectors on the n th processor in each group are denoted {g} n - 

Wavefunctions in real space are obtained by a parallel FFT, so real space points are distributed as for reciprocal 
space points, and this FFT is performed in parallel across all processors in each group. As efficient parallel FFT 
algorithms are available and a distributed sum is trivial these sections of the algorithm parallelise well. The integral 
must also be must be summed over pairs of k points, with each point in a different group, and it is this part of 
the algorithm that requires a large amount of waiting and communication, and is not as efficient. The algorithm for 
evaluation of Eq. (|16p takes the following form, with the current processor in group i. 



C*"(Pk) :=Transform[C*"(k)] 
«^(r) := FFT[C«(Pk)] 

for n' := l,AWds 
for m := 1, Nq 
for k' e {k}i 



III. IMPLEMENTATION AS A PARALLEL ALGORITHM 



[t] for n := 1, AT bands 
for k S {k} 4 



for S € Space-group 



If m = i then ug'(r) := FFT[C"'(k')] 
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Beast u£/(r), k', ww, e£, from group m to all groups 
Sum Ujpk(r) h.Vy u^! (r)e i ( k ' _pk ^ (r+a) over group i 
Evaluate (/(ej) - /(eg)) - eg!) 
Add contribution to r\ using weights 

-/(*&))*'(*£- 

end 

end 

end 

end 

end 

end 

Sum 77 over all groups 

The subprogram Transform: obtains the image of the wavefunction in reciprocal space, as defined by Eq. (|18p . 
It takes the following form, where the current processor is in the group associated with k, and holds vectors {g}i 

Transform: 

for m := 1, P{g} 

Beast {g,C"(k)} m from m th to all processors in group 
for g G {g} 2 

If p-ig g { g } m then Q(Pk) := C£_ lg (k) e -*( pk +s)- w 

end 

end 

return C£(Pk) 

It is useful to know the scaling of this algorithm with respect to the distribution of data across the processors, 
(Nq, P{g}), and the number of processors available N p = Nq x P{g}- The time taken for a reasonably large system 
is dominated by the time spent carrying out the FFT contained within the inner loop (or waiting for another group 
to broadcast the result of this FFT) and is given by 

i,| ac N s N* ands Nl^-P{ g }t FFT (P{g}) (21) 

where tpFT(M) is the time taken to carry out a parallel FFT on M processors. The speedup of the parallel FFT 
with increasing number of processors is complex and depends on the architecture of the parallel system [17j , but some 
general conclusions about how this effects the performance of the entire algorithm can be made. If a linear speedup 
of the FFT with respect to the available processors (P{g}) occured then the total run time would be independent 
of the distribution of the processors among k and g points. However, the actual speedup will be worse than linear 
due to the communication times, so for a given N p the best efficiency is achieved when the k points are distributed 
amongst as many processors as possible, or P{g} is as small as possible. 

IV. RESULTS 

Test calculations have been carried out for a H atom moving above the hep hollow site of a Cu(lll) surface. This 
system has been chosen for its simplicity (although we note that spin polarisation is needed to obtain the correct 
electronic structure at larger atom-surface separations) and because of its relevance to the chemicurrent experiment 
of Nienhaus et al The surface is modelled by a five-layer slab, with a vacuum gap equivalent to another five layers. 
A 2 x 2 in-plane super-cell is used - tests show that the deformation potential caused by the displacement of H atoms 
is well localised within this unit cell. A spin-polarised version of the PW91 functional is used for exchange-correlation 
effects ([H|, and references therein), a Troullier-Martins [l9| pseudopotential is used for Cu, and H is represented by 
a Coulomb potential. The plane-wave, pseudopotential code CASTEP is used to obtain the self-consistent potentials 
and Kohn-Sham states that are required for the calculation of the matrix elements in Eq. (fTo| . Calculations are 
performed with a plane- wave cut-off of 830 eV, 54 k-points are included in the full surface Brillouin zone, and a Fermi 
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T A T 


m + vi Vk= 


=k',T + Vk=k',l 


(K) (eV) 


(meV ps A 2 ) 


46 0.01 


1.578 


1.364 


231 0.05 


1.576 


1.361 


463 0.10 


1.570 


1.308 


926 0.20 


1.551 


1.217 


1389 0.30 


1.532 


1.166 


1852 0.40 


1.525 


1.146 


2315 0.50 


1.536 


1.151 



TABLE I: Variation of friction coefficient with temperature, A = 0.5 eV 

surface broadening of 0.25 eV is used. In order to test the convergence of the method we present results for the H 
atom 2.5 A above the surface, and with h perpendicular to the surface. 

To obtain the perturbative field, h.VV^, calculations are performed with H at x±h. Equation (fTU|) is then applied to 
the resulting self-consistent potentials to obtain h.W. The finite difference introduces two errors in the final result. 
First a quadratic error is introduced by the finite difference itself, and second any small errors in the Kohn-Sham 
potentials are magnified for small |h|. It follows that |h| must be carefully chosen to be small enough to minimise 
the first of these errors, but large enough to minimise the second. Tests for a number of |h| suggest that |h| = 0.02 
A results in a quadratic error w 0.1% and an error due to noisy Kohn-Sham potentials of w 1.0%. Quantifying the 
latter of these is not straightforward, hence a pessimistic estimate has been given. 

With these parameters Eq. (|16|) was evaluated using the algorithm described in section IIIII If the factor 

f/( e k) — /( e k')) ^'( e k — e k') m Eq. (HU) was less than 10~ 3 the contribution to r) was not calculated, increasing 
the efficiency of the calculation. Enough bands are included for the highest energy bands at each k point to be 
discarded. 

Two parameters remain, which describe the temperature of the system, and how the discretisation of k space is 
dealt with through the smearing, A. The temperature enters through the Fermi-Dirac occupation functions in Eq. 
(|16[) . and may be chosen to take any value. For the metallic system considered here a weak temperature dependence 
is expected, so the zero-temperature limit is the quantity of interest. We chose to use a 'squashed Fermi-Dirac' 
distribution, f(x), related to the function defined in Eq. (fT2]) by 

/x— ep 
5{x)dx. (22) 
-oo 

If At is the width associated with S(x) then the properties of the system are very close to those of a system described 
by a Fermi-Dirac distribution with a temperature given by Ax/(v27rfcB) ; where fee is the Boltzmann constant. 

Table [J shows results for r\ summed over both spins, calculated over a range of temperatures for A = 0.5 eV (this is 
discussed below) and for both the isolated atom (77^ see Eq. (fTo) ) and the coherent over-layer (?7k=k',T + 7 7k=k',j., 
see Eq. (flTj) K Results for an isolated atom are weakly dependent on temperature, and we take At = 0.1 eV to 
represent the low temperature limit. For the coherent over-layer the temperature dependence is stronger, which is 
understandable since a greater proportion of the available (band to band) transitions are expected to be of higher 
energy. 

The approximation that requires the most attention is the discretisation of k space and the reintroduction of a 
continuum through smearing. As described in section [TT] we derive a continuous self-energy function using the S(x) of 
Eq. (fT2|) with width A, and extract the linear behaviour of this function close to zero energy. For a given set of k 
points we must find a value of A large enough for the approximate self-energy to have converged to a linear function 
near zero energy. This critical A must be small enough (and so the density of the k point mesh must be high enough) 
to avoid the higher energy structure of the self-energy influencing the structure close to zero energy. 

Results were calculated as above, with At = 0.1 eV and A ranging from 0.05 to 1.0 eV. These are shown in Fig. 
[3] for both the isolated atom and coherent over-layer. For the isolated atom r\ shows good convergence by A = 0.6 
eV, and this value is low enough to conserve the structure of the self-energy. Convergence has not been achieved for 
the coherent over-layer, and the strong variation of the friction coefficient with A in this range suggests that the 54 
k-points are insufficient to achieve convergence for this system. This difference in the convergence behaviour is due 
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FIG. 3: Convergence of friction coefficient with the 'smearing parameter' A. Results are given for At = 0.1 eV (T = 463 K). 
The solid line shows the friction coefficient for an isolated atom above the surface, r\ = ry-p + r)i , whereas the dashed line shows 
the friction coefficient for an over-layer of atoms in coherent motion, ?7k=k',T + ^k=k',|- 



to far fewer transitions being available in the coherent over-layer system due to the k = k' condition, and so fewer 
transitions to approximate a continuous self-energy. 



V. CONCLUSION 



An ab initio method has been presented that allows the evaluation of the friction due to electron-hole pair creation 
experienced by an isolated molecule in motion near a metal surface. The approach described combines Kohn-Sham 
super-cell methods employing a plane-wave basis with a description of the electron-hole pair creation process via 
TDDFT. Results have been presented for a H atom above a Cu(lll) surface, and convergence of the calculations has 
been tested. We find a significant difference in the results for motion of an isolated atom and those for a coherent 
over-layer of atoms, both in the physical properties of the system and in their numerical calculation. Since the 
calculation is relatively expensive to perform for systems of interest an efficient parallel implementation of method 
has been given. 

We wish to thank M. Persson and S. Holloway for useful discussions. 
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